
####################
### ESS analysis ###
####################

rm(list = ls())
gc()

library(tidyverse)
library(stargazer)
library(data.table)
library(lfe)
library(BBmisc)
library(lme4)

### Set WD
setwd("~/Dropbox (Princeton)/Research/*Spatial_Inequality/paper_SMD_inequality/JoP_final/replication")  

### Multilevel data for ESS 9 
load("data/ess9_ml.RData")

### Prepare
ess9_ml[ hinctnta == "J - 1st decile" , hh_inc := 1 ]
ess9_ml[ hinctnta == "R - 2nd decile" , hh_inc := 2 ]
ess9_ml[ hinctnta == "C - 3rd decile" , hh_inc := 3 ]
ess9_ml[ hinctnta == "M - 4th decile" , hh_inc := 4 ]
ess9_ml[ hinctnta == "F - 5th decile" , hh_inc := 5 ]
ess9_ml[ hinctnta == "S - 6th decile" , hh_inc := 6 ]
ess9_ml[ hinctnta == "K - 7th decile" , hh_inc := 7 ]
ess9_ml[ hinctnta == "P - 8th decile" , hh_inc := 8 ]
ess9_ml[ hinctnta == "D - 9th decile" , hh_inc := 9 ]
ess9_ml[ hinctnta == "H - 10th decile" , hh_inc := 10 ]

ess9_ml[ gndr == "Male" , female := 0] 
ess9_ml[ gndr == "Female" , female := 1] 

ess9_ml[ , married_fin := ifelse(maritalb %in% c("Legally married","In a legally registered civil union"), 1, 0)]
ess9_ml[ is.na(maritalb), married_fin := NA ]

ess9_ml[ , has_children := ifelse(chldhhe %in% c("Yes"), 1, 0)]
ess9_ml[ is.na(chldhhe), has_children := NA ]

ess9_ml[ , unemployed := ifelse(uempla == "Marked" , 1, 0 )]
ess9_ml[ , retired := ifelse(rtrd == "Marked" , 1, 0 )]
ess9_ml[ , employed := ifelse(pdwrk == "Marked" , 1, 0 )]
ess9_ml[ , student := ifelse(edctn == "Marked" , 1, 0 )]

ess9_ml[ , type_living_area := domicil ]
ess9_ml[ domicil %in% c("No answer","Don't know","Refusal") , type_living_area := NA ]

ess9_ml[ , type_living_area2 := type_living_area ]
ess9_ml[ type_living_area %in% c("Country village","Farm or home in countryside") , type_living_area2 := "village/farm" ]

ess9_ml[ , urban := ifelse(type_living_area %in% c("A big city"), 1, 0) ]
ess9_ml[ is.na(type_living_area), urban := NA ]


ess9_ml[ lrscale %in% seq(1,9,1) , lrscale_v2 := lrscale ]
ess9_ml[ lrscale %in% c("Left") , lrscale_v2 := "0" ]
ess9_ml[ lrscale %in% c("Right") , lrscale_v2 := "10" ]
ess9_ml[ , lrscale_v2 := as.numeric(as.character(lrscale_v2)) ]

ess9_ml[ gincdif == "Agree strongly" , gvt_reduce_ineq := 2]
ess9_ml[ gincdif == "Agree" , gvt_reduce_ineq := 1]
ess9_ml[ gincdif == "Neither agree nor disagree" , gvt_reduce_ineq := 0]
ess9_ml[ gincdif == "Disagree" , gvt_reduce_ineq := -1]
ess9_ml[ gincdif == "Disagree strongly" , gvt_reduce_ineq := -2]

ess9_ml[ , gvt_reduce_ineq_norm := BBmisc::normalize(gvt_reduce_ineq, method="range", range=c(0,1)) ]

ess9_ml[ , foreign_born := ifelse(brncntr %in% c("No"), 1, 0)]
ess9_ml[ is.na(brncntr), foreign_born := NA ]

ess9_ml[ !is.na(ess9_reg) , count := .N, by=c("cntry","ess9_reg")]

ess9_ml[ eisced %in% c("ES-ISCED I , less than lower secondary"), eisced_v2 := "ES-ISCED I" ]
ess9_ml[ eisced %in% c("ES-ISCED II, lower secondary"), eisced_v2 := "ES-ISCED II" ]
ess9_ml[ eisced %in% c("ES-ISCED IIIa, upper tier upper secondary","ES-ISCED IIIb, lower tier upper secondary"), eisced_v2 := "ES-ISCED III" ]
ess9_ml[ eisced %in% c("ES-ISCED IV, advanced vocational, sub-degree"), eisced_v2 := "ES-ISCED IV" ]
ess9_ml[ eisced %in% c("ES-ISCED V1, lower tertiary education, BA level","ES-ISCED V2, higher tertiary education, >= MA level"), eisced_v2 := "ES-ISCED V" ]
ess9_ml[ eisced %in% c("Other"), eisced_v2 := "ES-other" ]
ess9_ml[ , eisced_v2 := factor(eisced_v2, levels=c("ES-other","ES-ISCED I","ES-ISCED II","ES-ISCED III","ES-ISCED IV","ES-ISCED V")) ]

ess9_ml[ , hh_inc_fac := as.factor(hh_inc) ]
ess9_ml[ , hh_inc_fac := relevel(hh_inc_fac,"1")]

### Remove countries with missing data and less than 3 regions
ess9_ml_reg_ineq <- ess9_ml[ !is.na(ess9_reg) & count >= 50 & ! cntry %in% c("AL","ME","RS","XK","LU","CY","IS") & !is.na(hh_inc) ]

### Inequality measures 
ess9_ml_reg_ineq[ !is.na(ess9_reg), reg_gini_hh_inc_wght := lapply(.SD, function(x) wINEQ::Gini(x, W=dweight) ) , by=c("cntry","ess9_reg"), .SDcols="hh_inc" ]

### Prep
ess9_ml_reg_ineq[ , N_cnty := .N, by=c("cntry")]

d1 <- unique(ess9_ml_reg_ineq[ ,.(cntry,ess9_reg)])
d1.1 <- d1[ , .N, by=c("cntry")]

d2 <- unique(ess9_ml_reg_ineq[ !is.na(n2_gdp_eurhab_2019) ,.(cntry,ess9_reg)])
d2[,.N, by=cntry]

d3 <- unique(ess9_ml_reg_ineq[ ,.(cntry,regunit,N_cnty)])
d3.1 <- merge(d3,d1.1,by="cntry")

d3.2 <- d3.1[,.(cntry,N,regunit,N_cnty)]

d3.2[ , regunit := gsub("NUTS level ","",regunit) ]

names(d3.2) <- c("Country","Number of regions","NUTS level","Total number of respondents")

### SI Table C.1
print(xtable::xtable(d3.2, include.colnames=T, caption="Summary statistics"), include.rownames=F)





### ------------
### Models 
### ------------

mod1.1 <- felm(gvt_reduce_ineq_norm ~ reg_gini_hh_inc_wght  | cntry | 0 | ess9_reg, ess9_ml_reg_ineq, weights = ess9_ml_reg_ineq$dweight)

mod1.2 <- felm(gvt_reduce_ineq_norm ~ reg_gini_hh_inc_wght + female + agea + foreign_born + eisced_v2 + employed + unemployed + student + retired + hh_inc_fac + 
                  married_fin + has_children | cntry | 0 | ess9_reg, ess9_ml_reg_ineq, weights = ess9_ml_reg_ineq$dweight)

mod1.3 <- felm(gvt_reduce_ineq_norm ~ reg_gini_hh_inc_wght + female + agea + foreign_born + eisced_v2 + employed + unemployed + student + retired + hh_inc_fac + 
                 married_fin + has_children + urban + lrscale_v2 | cntry | 0 | ess9_reg, ess9_ml_reg_ineq, weights = ess9_ml_reg_ineq$dweight)

mod1.4 <- felm(gvt_reduce_ineq_norm ~ reg_gini_hh_inc_wght + female + agea + foreign_born + eisced_v2 + employed + unemployed + student + retired + hh_inc_fac + 
                 married_fin + has_children + urban + lrscale_v2 + 
                 n2_tpopsz_2020 + n2_pode_2019 + n2_mio_eur_2019 + n2_unraall_2020 | 0 | 0 | ess9_reg, ess9_ml_reg_ineq, weights = ess9_ml_reg_ineq$dweight)

mod1.5 <- lmer(gvt_reduce_ineq_norm ~ reg_gini_hh_inc_wght + female + agea + foreign_born + eisced_v2 + employed + unemployed + student + retired + hh_inc_fac + 
                 married_fin + has_children + urban + lrscale_v2 + 
                 n2_tpopsz_2020 + n2_pode_2019 + n2_mio_eur_2019 + n2_unraall_2020 + (1|cntry) + (1|ess9_reg), ess9_ml_reg_ineq, weights = ess9_ml_reg_ineq$dweight)

mod1.6 <- felm(gvt_reduce_ineq_norm ~ reg_gini_hh_inc_wght + female + agea + foreign_born + eisced_v2 + employed + unemployed + student + retired + hh_inc_fac + 
                 married_fin + has_children + urban + lrscale_v2 | cntry | 0 | ess9_reg, ess9_ml_reg_ineq[ regunit == "NUTS level 3" ], weights = ess9_ml_reg_ineq[ regunit == "NUTS level 3" ]$dweight)


covar.lab <- c("Regional inequality","Female","Age","Foreign born","ES-ISCED I","ES-ISCED II","ES-ISCED III","ES-ISCED IV","ES-ISCED V",
               "Employed","Unemployed","Student","Retired","HH income decile 2","HH income decile 3","HH income decile 4",
               "HH income decile 5","HH income decile 6","HH income decile 7","HH income decile 8","HH income decile 9","HH income decile 10",
               "Married","Has children","Urban residence","Left-right ideology","NUTS2 Population","NUTS2 Population density","NUTS2 GDP (mn. EUR)","NUTS2 Unemployment rate")


### SI Table C.2
stargazer(mod1.1, mod1.2, mod1.3, mod1.4, mod1.5, mod1.6, 
          digits.extra=0, digits=2,  
          align=T, no.space=T, omit.stat=c("ser","F"), type="text",
          dep.var.labels = c("Should the government reduce income inequality"),
          covariate.labels = covar.lab, omit = "Constant",
          title = "Effect of Regional Inequality on Support for Redistribution",
          table.placement = "h!", 
          add.lines = list(c("Mean DV",
                             round(mean(mod1.1$fitted.values),2),
                             round(mean(mod1.2$fitted.values),2),
                             round(mean(mod1.3$fitted.values),2),
                             round(mean(mod1.4$fitted.values),2),
                             round(mean(fitted(mod1.5)),2),
                             round(mean(mod1.6$fitted.values),2)),
                           c("Country Fixed Effects","$\\checkmark$","$\\checkmark$","$\\checkmark$","-","-","$\\checkmark$"),
                           c("Country Random Effects","-","-","-","-","$\\checkmark$","-")))


